source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
View(slaves)
library(poweRlaw)
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
p = rpldis(10^6, 1, 2.5 )
bins = slaves$Slaves_low + 1
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
slaves$Slaves_low = Slaves_low + 1
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
slaves$Slaves_low = slaves$Slaves_low + 1
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
bins = c(slaves$Slaves_low, 10^9)
p = rpldis(10^6, 1, 2.5 )
h = hist(p, breaks = bins, plot = F)
View(slaves)
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
View(slaves)
h = hist(p, breaks = bins, plot = F)
h = hist(p, breaks = bins)
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
View(slaves)
plot(slaves$Slaves_low, slaves)
plot(slaves$Slaves_low, slaves$share)
plot(slaves$Slaves_low, slaves$share, log = "xy")
plot(slaves$Slaves_low, slaves$share, log = "xy")
points(slaves$Slaves_low, slaves$share, log = "xy")
points(slaves$Slaves_high, slaves$share, log = "xy")
library(zoo)
mids = rollmean(bins)
bins = c(slaves$Slaves_low, 10^9)
bins = c(slaves$Slaves_low, 10^4)
mids = rollmean(bins)
p = rpldis(10^6, 1, 2.5 )
bins = c(slaves$Slaves_low, 10^4)
mids = rollmean(bins)
mids = rollmean(bins,2)
p = rpldis(10^6, 1, 2.5 )
h = hist(p, breaks = bins, plot = F)
plot(mids, slaves$share)
plot(mids, slaves$share, log ="xy")
points(mids, h$density)
points(mids, h$counts/sum(h$counts))
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
mean(p)
bins = c(slaves$Slaves_low, Inf
b_min = min(bins)
n = sum(bin_sub$counts)
h = bin_sub$counts
b_l = bin_sub$lower
b_u = bin_sub$upper
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
bins = c(slaves$Slaves_low, max(p))
p = rpldis(10^6, 1, 2.7 )
mids = rollmean(bins,2)
h = hist(p, breaks = bins, plot = F)
plot(mids, slaves$share, log ="xy")
points(mids, h$counts/sum(h$counts), col = "red")
bins = c(slaves$Slaves_low, Inf)
b_min = min(bins)
n = sum(slaves$Number_owners)
b_u = Slaves_high
b_u = slaves$Slaves_high
b_min = min(bins)
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = slaves$Slaves_low
b_u = slaves$Slaves_high
View(slaves)
b_min = 1
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = slaves$Slaves_low
b_u = slaves$Slaves_high
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
n/10^6
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
tail(b_u,1) = Inf
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
tail(b_u,1) = Inf
b_u
tail(b_u,1)
tail(b_u,1) = Inf
tail(b_u,1) = Inf
b_u[length(b_u)] = Inf
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
b_min = 1
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = slaves$Slaves_low
View(slaves)
View(slaves)
slaves = fread("Atibellum_slave_dist.csv")
View(slaves)
b_min = 1
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = slaves$Slaves_low
b_u = slaves$Slaves_high
b_l[1] = 1
b_u[length(b_u)] = Inf
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
b_min = 1
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = slaves$Slaves_low
b_u = slaves$Slaves_high
b_u[1] = 1
b_u[length(b_u)] = Inf
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
bins = c(slaves$Slaves_low, Inf
b_min = 1
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = slaves$Slaves_low
b_u = slaves$Slaves_high
b_u[1] = 1
b_u[length(b_u)] = Inf
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
bins = c(slaves$Slaves_low, max(p))
p = rpldis(10^6, 1, 2.7 )
mids = rollmean(bins,2)
h = hist(p, breaks = bins, plot = F)
plot(mids, slaves$share, log ="xy")
points(mids, h$counts/sum(h$counts), col = "red")
bins = c(slaves$Slaves_low, Inf)
b_min = bins[1]
library(here)
library(ineq)
library(data.table)
library(poweRlaw)
library(zoo)
dir = here()
setwd(dir)
slaves = fread("Atibellum_slave_dist.csv")
slaves$Slaves_low = slaves$Slaves_low + 1
slaves$Slaves_high = slaves$Slaves_high + 1
slaves$share = slaves$Number_owners/sum(slaves$Number_owners)
bins = c(slaves$Slaves_low, Inf)
b_min = bins[1]
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = bins[-length(bins)]
b_u = bins[-1]
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
bins = c(slaves$Slaves_low, max(p))
p = rpldis(10^6, 1, 2.72 )
mids = rollmean(bins,2)
h = hist(p, breaks = bins, plot = F)
plot(mids, slaves$share, log ="xy")
points(mids, h$counts/sum(h$counts), col = "red")
mean(p)
p = rpldis(10^6, 1, fit$par )
h$density
# power law random sample
#####################################################
library(doSNOW)
library(tcltk2)
library(snow)
library(foreach)
cl <- makeCluster(n.core, type="SOCK")
cl <- makeCluster(8, type="SOCK")
registerDoSNOW(cl)
clusterSetupRNG (cl, type = "RNGstream")
beginTime <- proc.time()
pb <- txtProgressBar(max = n.iteration, style = 3)
n.iterarion = 10
beginTime <- proc.time()
pb <- txtProgressBar(max = n.iteration, style = 3)
n.iteration = 10
# Progress Bar
beginTime <- proc.time()
pb <- txtProgressBar(max = n.iteration, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress=progress)
# 6. Run model for each set of parameters  ----------------------------------------------------
model =  foreach (iteration = 1:n.iteration, .options.snow=opts, .combine=rbind,
.packages=c('ineq', 'poweRlaw')) %dopar%  {
p = rpldis(10^6, 1, fit$par )
bins = c(slaves$Slaves_low, max(p))
h = hist(p, breaks = bins, plot = F)
h$density
}
View(model)
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
library(hmod)
p = rpld(10^6, 1, fit$par, discrete_max = 10^5, ordered = T)
p = rpldis(10^6, 1, fit$par )
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
View(slaves)
View(model)
conf = function(x){
p = quantile(x, probs = c(0.025, 0.5, 0.975))
}
test = apply(model, 2, conf)
View(test)
test = t( apply(model, 2, conf) )
conf_interval = t( apply(model, 2, conf) )
result = data.table(slaves$Slaves_low, Slaves_high, conf_interval)
result = data.table(slaves$Slaves_low, slaves$Slaves_high, conf_interval)
View(result)
conf = function(x){
p = quantile(x, probs = c(0.05, 0.5, 0.95))
}
conf_interval = t( apply(model, 2, conf) )
result = data.table(slaves$Slaves_low, slaves$Slaves_high, conf_interval)
result = data.table(slaves$Slaves_low, slaves$Slaves_high, conf_interval)
View(result)
View(model)
View(slaves)
View(result)
bins
slaves$density = slaves$Number_owners / diff(bins)
View(slaves)
p = rpld(10^6, 1, fit$par, discrete_max = 10^5, xmax = 10^5, ordered = T)
bins = c(slaves$Slaves_low, 10^5)
stopCluster(cl)
stopCluster(cl)
slaves$density = slaves$Number_owners/diff(bins)
View(slaves)
sum(slaves$density)
slaves$density = slaves$Number_owners / diff(bins) / sum(slaves$Number_owners)
sum(slaves$density)
sum(slaves$density*diff(bins))
result = data.table(slaves, conf_interval)
View(result)
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Soltow/slave_dist_fit.R')
View(result)
plot(result$Slaves_low, result$density)
plot(result$Slaves_low, result$density, log ="xy")
points(result$Slaves_low, result$`5%`)
points(result$Slaves_low, result$`5%`, col = "red")
points(result$Slaves_low, result$`95%`, col = "red")
points(result$Slaves_low, result$`50%`, col = "red")
fwrite(results, "slaves_power_law_fit.csv")
fwrite(result, "slaves_power_law_fit.csv")
library(here)
library(ineq)
library(data.table)
library(poweRlaw)
library(zoo)
dir = here()
setwd(dir)
# https://www.measuringworth.com/slavery.php
slaves = fread("Atibellum_slave_dist.csv")
# power law fit
###########################################################
bins = c(slaves$Slaves_low, Inf) + 1
b_min = bins[1]
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = bins[-length(bins)]
b_u = bins[-1]
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
slaves = fread("Antebellum_slave_dist.csv")
library(here)
library(ineq)
library(data.table)
library(poweRlaw)
library(zoo)
dir = here()
setwd(dir)
# data from https://www.measuringworth.com/slavery.php
slaves = fread("Antebellum_slave_dist.csv")
# power law fit
###########################################################
bins = c(slaves$Slaves_low, Inf) + 1
b_min = bins[1]
n = sum(slaves$Number_owners)
h = slaves$Number_owners
b_l = bins[-length(bins)]
b_u = bins[-1]
log_lik = function(alpha){
-(  n*(alpha - 1) * log(b_min) + sum( h*log( b_l^(1-alpha) - b_u^(1-alpha) ) ) )
}
fit = nlminb(3, log_lik, lower = 1.5, upper = 9)
# power law random sample
#####################################################
library(doSNOW)
library(tcltk2)
library(snow)
library(foreach)
library(hmod)
n.iteration = 1000
cl <- makeCluster(8, type="SOCK")
registerDoSNOW(cl)
clusterSetupRNG (cl, type = "RNGstream")
beginTime <- proc.time()
pb <- txtProgressBar(max = n.iteration, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress=progress)
model =  foreach (iteration = 1:n.iteration, .options.snow=opts, .combine=rbind,  .packages=c('ineq', 'poweRlaw', 'hmod')) %dopar%  {
p = rpld(10^6, 1, fit$par, discrete_max = 10^5, xmax = 10^5, ordered = T)
# shift p distribution to start at 0
p = p - 1
bins = c(slaves$Slaves_low, max(p))
h = hist(p, breaks = bins, plot = F)
h$density
}
stopCluster(cl)
# model confidence intervals
#######################################################
conf = function(x){
p = quantile(x, probs = c(0.05, 0.5, 0.95))
}
conf_interval = t( apply(model, 2, conf) )
# model confidence intervals
#######################################################
conf = function(x){
p = quantile(x, probs = c(0.05, 0.5, 0.95))
}
conf_interval = t( apply(model, 2, conf) )
# slave density
bins = c(slaves$Slaves_low, 10^5)
slaves$density = slaves$Number_owners / diff(bins) / sum(slaves$Number_owners)
# write results
result = data.table(slaves, conf_interval)
fwrite(result, "slaves_power_law_fit.csv")
View(result)
result = data.table(slaves, model = conf_interval)
View(result)
result = data.table(slaves, power_law = conf_interval)
View(result)
conf = function(x){
p = quantile(x, probs = c(0.05, 0.5, 0.95))
}
conf_interval = t( apply(model, 2, conf) )
# slave density
bins = c(slaves$Slaves_low, 10^5)
slaves$density = slaves$Number_owners / diff(bins) / sum(slaves$Number_owners)
# write results
result = data.table(slaves, power_law = conf_interval)
fwrite(result, "slaves_power_law_fit.csv")
View(conf_interval)
View(model)
View(result)
